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1. Introduction 

Following the ideas of A.N. Tykhonov about ill-posed problems inversion, we extend the 
iterative regularization method by including alternatively between two iterations least 
square gradient looking terms, and with a modification of the parametrization of the 
regularizing term. The merits of this approach are a fast convergence to a satisfying 
solution, stable and unsensitive to the noise fluctuations, and a complete independence 
with respect to the initial guess. With this method, we unfold one dimensional data 
set of measurements with a Bonner Sphere system to obtain a wide neutron spectrum 
estimation. Two experimental cases are presented : a complex thermal neutron wide 
spectrum and a single monoernegetic 565 keV neutron spectrum. The evolution of 
the criteria parameter is discussed and also the systematics and statistics error bars 
calculations. 

2. Theoretical approach 

Once we have a set of 11 measurements, by spline interpolation, we obtain a vector 
y of dimension 64 {y G M^^"^); this is justified by the smoothness of the data set 
and by the convolution that erodes most of the details of the initial spectrum. The 
measurement function y G C°° and is very smooth. The spline interpolation produces a 
64 bin vector, with continous derivatives at each of the 11 knots and local second degree 
polynomial interpolation. These 11 measurements corresponds to the responses to a 
neutron spectrum of an Helium^ neutron detector with various polyethylene spheres, 
that are equivalent to different filters with response matrix A. With x G R+^^ as the 
unknown neutron spectrum, y as the measurements, the convolution equation of the 
system stands as : 



As the problem is ill-posed, neither the direct inversion, nor the least square solution are 
satisfying solutions, because they are irregular and highly sensitive to the fiuctuations 
of the measurement's uncertainties. That means that continuous variations of y does 
not lead to continuous variations in x. 



y — Ax 
X G 



(1) 



y G A G X 



Direct inversion : 



estimator : x — A~ y 
Least square solution : 

criterion : Q = \Ax — y 



(2) 



estimator : x = [A^A] ^A^y 
(A*; Transposed matrix) 



(3) 
(4) 
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2.1. Iterative Tykhonov approach 

A better solution can be achieved through the scheme of the basic Tykhonov method 
[SI m El El IHl Ej , which is more regular and less sensitive to noisy fluctuations. 
A regularization term is added to the least square criterion through a Lagrange 
multiplicator method. This term ensures the continuity of the solution and its closedness 
to an initial guess Xq, which can be a "not so bad" first solution. 

Qx = \Ax - + X\x - xol"^ Ae[0,+oo] (5) 
x=[A'A + XI]-^ [A'y + Xxo] (6) 

The demonstration of the calculation of the estimator (equation 6) can be found as 
follows P : 

Assume x to be the global minimum of the criterion Qx- 

Vx G M+'^^ : Qx,, < Qx,, (7) 
This is true, moreover in the case of small variations around x : 

x = x + ah he a G M 

Q\,x < Q\,a,h 

Q\,a,h — Qx,x > 

■x^ \A{x + ah) — + X\x + ah — xqI'^ — \Ax — — A|2; — a;oP > 
Due to the definition of the scalar product : 

^ a*{AhY[Ax -y] + a[Ax - yfAh + \a\'^\Ah\'^ + a*Xh^[x - Xq] + 
aX[x - xoYh + X\a\^\h\'^ > 

{a*h^){A^[Ax -y] + X[x - xo]) + a{[AS: - yfA + A[£ - xoY)h > 
That is independent of /i* and a ; then it implies : 
=^ (A^Ax - A^y + Xx - Xxq) = 
^ [A^A + XI]x - [A^y + Xxo] = 
^x = [A' A + XI]-^[A'y + Xxo] 

The uniqueness and existence of a bounded solution is also proved in Bertero P3I21- The 
Tykhonov iterative approach is achieved by successive refinement of the solution with 
several iterations, beginning with an initial guess Xq, the solution xi satisfies better the 
criterion of Tykhonov than xq, and so on : 

Qx,n = \AXn - yl"^ + X\Xn - Xn-l\'^ A G [0, +Oo] (8) 

Xn = [A'A + XI]~'[A'y + Xxn-i] (9) 

Q\,l > Q\,2 > ••• > Q\,n > Q\,n+1 > ••• 

2.2. Modification of the running parameter X 

In order to help the computer calculation, a slight modification is done with the 
modification of the running interval of the parameter A. We also decided that the 
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two terms of the Qx^n criterion have to be of equal importance : 

Ag[0,+oo]i — ^ AG [0,1] (10) 

Qx,n = {l-X)\AXn-y\^ + \\Xn-Xn-l\^ A G [0, l] (11) 

Xn = [(1 - X)A'A + A/]-i[(l - X)A'y + Ax„„i] (12) 
Let's have Xq as an initial guess, and lets compute Xi. The case A = gives : 

Qo^i = \Axi - y\^ (13) 

=^ xo,i = [A^A]^^A^y (least square solution) (14) 

^ < go,i (15) 

On the other side, the case A = 1 gives : Qi.i = |xi — xop with the solution xi = xq then 
Qi i = then A = 1 is a minimorum and all the iterations stand for an identity operator, 
then Wn G N* Xn = xq independently of xq ! So this method gives always a wrong 
solution, because all initial guesses are available solutions without any restriction. We 
fixed this problem by modifying the criteria in the following way : 

Q^^^= i^-mxn-y\' + XK-x...,\^ ^^^^^^^ ^^g^ 



:i-A)A 

Qx,n = \ YZTx — ^ 



the estimator is unchanged : 

x„ = [(1 - X)A'A + A/]-i[(l - X)A'y + Aa;„_i] (18) 
but Qo^n = +00 and = +oo 

Also we can note : 3A g]0, 1[ , Q\^n < +oo -v^ 3 an actual minimorum of Qx^n and the 
false minimorum Qi^n disappears! 

2.3. Possible generalization to derivatives of order m 

Basically, the Tykhonov method can be generalized to derivatives of order m of Xn 

instead of using the simple squared modulus [TUI. These derivatives D^"^' 

the derivatives of order m of the spectrum Xn{E) regarding the energy binning E : 

D^X, = ^Xn{E) 

This generalization to greater orders of derivatives is given by : 

lAxn-yl"^ , , 

Qx,n,m = ' \ ^' + ^Y^A ^ ^ 



Xr 



[(1 - X)A'A + AD('")*D(™)]-1[(1 - X)A'y + XD^'^'^'D^'^^Xn-i] (20) 



But we have found that the convergence to a stable solution is rather better and 
faster with a simple function without the m^^ derivative, so we decided to restrict our 
calculation only to order 0. 
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2.4- Hybridation with the least square method 

Three main concerns guided the development of this unfolding iterative method : a fast 
converging method, a better solution than the classical Tykhonov method or than the 
least square method and the last one, a complete independence of the solution regards to 
the initial guess. The iterative Tykhonov approach ensures the last concern but not the 
first two. So we decided to introduce some accelerating (focusing) algorithm between 
two classical Tykhonov iterations. Let suppose that Xn-i and x„ are two successive 
solutions obtained with this Tykhonov algorithm, we create some kind of derivative 
term regarding the iteration parameter n: 



lim — ^ = 0^3xE M+^^; lim Xn = x 

on n^+oo 



dXn/dn = Xn- Xn-l 

dx 

lim — ^ = ^ 3a 

n^+oo On 

Then we minimize Q^ n+i 

doc 

Qf,,n+i = \Axn+i - yf with Xn+1 = x„ + Z^"^; G M 

^Q^,,n+l = \A{Xn+fl-^)-y\^ (21) 

So after every two Tykhonov unfoldings, we compute one solution closed to the least 
square solution. 

At the end of the convergence, the reconstructed spectrum = Axn is closer to the 
measured spectrum y than with use of the classical Tykhonov iterative scheme, which 
presents the compromise of an enhancement of the stability of the solution toward the 
noise fluctuations while the price to pay is a loss of the accuracy ^ and an excursion 
from the measured spectrum y. By the way, this approach converges faster than a 
classical iterative Tykhonov method. This hybridation algorithm acts as a focalisation 
method with respect to the requested measurement y, and also as a great accelerator of 
convergence. 

2.5. Iteration schema 

The iterative schema can be described in the following way. Let's have xq as any initial 
guess. 

iteration 1 : 

\Axi-y\'^ -xqP 

A + (1-A) ^^^^ 

xi = [(1 - X)A'A + XI]-^[{1 - X)A'y + Axq] (23) 

iteration 2 : 

_ \Ax2-y\^ \x2-xi\^ 

^^■^ A + (1-A) ^^^^ 

X2 = [(1 - X)A'A + XI]-^[il - X)A'y + Xx,] (25) 

iteration 3: 
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dn 



Q^^,= \A{x2 + f^^)-y\' (26) 



iteration 4 



0x2 

X3 = X2+ (27) 



\Ax^-y\^ IX4-X3I' 
QxA = ^ + (28) 

X4 = [(1 - A)A*A + A/]-^[(l - \)A'y + Xxs] (29) 



iteration n : n = 1 modulo 3 

\AXn y\ \Xn /■nn\ 
VA,n — \ 1 - \ ^ ^ 

x„ = [(1 - X)A'A + A/]-i[(l - X)A'y + Xx^-i] (31) 
iteration n+1 : n + 1 = 2 modulo 3 

VA,n+l - -j^ \ 

Xn+i = [(1 - X)A'A + A/]-i[(l - X)A'y + Ax„] (33) 
iteration n+2 : n + 2 = modulo 3 

QM,n+2 = |^(a;„+i + 1^—^) - y? (34) 



We have seen with the help of a computer that in all cases, this iterative schema 
converges to one unique solution, but we have not build the mathematical proof of this 
convergency! Practically, on a DEC-Compaq Alpha Workstation working with double 
precision 64 bits numbers, the resulting estimators are always stable after 300 iterations 
( 1 minute of computing time) and oftenly around 80 iterations, it is already stable. 
That is why we decided to put a cut on the maximum number of iterations at 300. We 
also decided to restrict the binning (in energy and interpolated sphere) over 64 bins, 
because working with 32 or 128 bins gives either a too bad energy resolution or a too 
slow computation, thus 64 seems to be an optimum. 



2. 6. Positivity of the solution 

Another important aspect of the unfolding is treated without any subtility: positivity 
|10j . Of course each bin of the spectrum must have a positive content, on the other case 
the spectrum would be unphysical. The positivity is introduced abruptly at each step 
of the iteration with the following non-linearity: 

Xn,i^ — ^ ^^^[O,^^,^ (36) 
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2. 7. Calculation of the statistics and systematics fluctuations 

The statistics fluctuations are estimated with a Monte-Carlo method. For one set of 11 
measurements with the Bonner Sphere system, 11 uncertainties of these measurements 
are calculated taking into account the number of observed counts and the duration 
of the measures. The propagation through the unfolding process of the statistics 
fluctuations of the measurements is computed with 100 sets of 64 data y, which are 
interpolated with the spline method from these 1 1 knots. The random generation of the 
11 knots is performed with a gaussian law centered onto the 11 actual measurements 
and with a standard deviation corresponding to the error bars directly measured with 
the instrument. The unfolding of the 100 data sets y gives 100 x estimators of the 
neutron spectrum. The fluctuations of the 100 spectra x correspond to the propagation 
of the statistics fluctuations through the unfolding process. 

j, statistics = ]y ^ = i^') 

Also, to the 100 x spectra correspond 100 reconstructed measurements y — Ax. The 
reconstructed measurements y are reinjected into the unfolding process giving new 
estimators of the spectra Xnew and twice reconstructed measurements ijnew The average 
of the absolute difference of the old and of the new spectra gives the systematics errors 
of the unfolding spectra : 

<systemat^cs = '"^^'^ "^"'^^^ ^ = (38) 

These two calculations of uncertainties are performed with 100 samples of Monte-Carlo 
generated data sets; this figure seems low due to the signal to noise ratio of 10 %, 
nevertheless the precision is sufficient for our purpose and moreover it is highly time 
consuming and we cannot afford greater computation. Also some tricks are used to avoid 
spending too much time recalculating 100 times almost the same things. Of course the 
A and /j, parameters minimizations are very expensive on the CPU time balance, so the 
evolution with the iterations of the A„ (for example) is memorized during the first pass, 
while for the next passes the minimization of the A parameter is performed at the n^^ 
iteration in the neighboorhood of the A„ first evaluation. Another trick is used, oftenly 
after 50 or 80 iterations the x„ and y„ figures are stables due to the precision of the 
computation, so each successive iteration appears as an Identity operator, then at this 
moment the iterative process is stopped. 



2.8. The initial guess 

This kind of unfolding method shows several advantages. The most important one is 
the independance of the final solution regarding the initial guess. Thus any distribution 
can be used as an initial guess. We have tested several different initial guesses : random 
numbers, flat functions Xq = 1 and the intercorrelation product between the response y 
and the discrete response functions of matrix A. 

xqj =< yi, Aij6j > 6j : dirac function (39) 
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The strong point of this method is that using all these different initial guesses the 
iterative approach converges to exactly the same result. For a reason of fast and easy 
computation, we decided to use a flat initial guess : xq — 1- 

2.9. Remarks about the minimization algorithms 

To circonvene the problem of the minimization with A of the Q{\) criteria, we used 
different technics in order to maintain a high precision on the real A and a reasonable 
time of calculation for each iteration. The minimization of the criteria is performed in 
different ways to optimize this calculation time. Firstly, a multiple step scanning and 
secondly, a dichotomy algorithm are performed. The first solution of the first set y is 
calculated carefully, while the other solutions computed to estimate the systematic and 
the statistic uncertainties, does begin their minimization at each step, with initial figures 
in the neighborhood of the first A„ sequence calculated previously. Of course, at this 
point the focusing on the right A is performed with scanning and dichotomy algorithms 
but with less steps, that avoid to spoil too much computation time. 

3. Experimental approach 

To illustrate this unfolding technics, we presents two examples of unfolding: the first 
one is a thermal neutron spectrum, the second one is a 565 keV mono-energetic neutron 
spectrum. 

3.1. Unfolding a thermal neutron spectrum from the Sigma installation 

This thermal neutron spectrum is produced near the Sigma facility at (IRSN- 
Cadarache), it consists of an assembling of several AmBe neutron sources surrounded 
by some graphite shielding. The expected spectrum is an heavy thermal component 
with some residual peaks between 1 and 3 MeV. Figures l.a and l.b shows precisely 
what was expected : a huge thermal neutron part and a small high energy component 
that goes from 100 keV up to 10 MeV. If the lower limit is physically explainable due 
to the tail of the decelarated neutrons from the high energy peaks, the upper bound 
around 10 MeV is not realistic and it is just due to the intrinsic resolution of the Bonner 
Sphere method itself. While figure l.a shows the statistics error bars, figure l.b presents 
the systematics unfolding error bars; in both case they are extremely tiny : the global 
statistics uncertainties is 0.11% while the global systematics errors of the unfolding is 
3.6%. These very small uncertainties are mainly due to the very high number of counts 
registered for most of the sphere, this can be seen on figure l.c, which represents the 
measured counts on the various sphere, (plain line), while the dot line corresponds to 
the recalculated spectrum after unfolding. These two shapes are in perfect agreement 
together. The last figures, shows one of the very specific parameter of the evolution 
of the unfolding process : 1 — \ vs. iteration number. This shows that after only 
20 iterations this parameter is already very low (10~^^), that means that the iterative 
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operator is very closed to the Identity operator, and the current solution is very closed to 
the final solution. This shape is very easy to unfold, this is why after only 78 iterations, 
the unfolding process is over. The main aspect of this unfolding is its ability to reveal 
some realistic details of the shape here a double structured shape with confortable 
uncertainties. 

3.2. Unfolding a monoenergetic neutron spectrum at 565 keV 

The Bonner Sphere system was operated in the beam of a Van de Graaf accelerator 
producing 565 keV mono-energetic neutrons at CEA (Bruyeres le Chatel). The expected 
spectrum is a single energy peak enlarged by the energy resolution of the device 
eventually surrounded by some low energy diffused neutrons. The background of the 
low energy diffused neutrons is substracted from the measurements thanks to other 
measures undertook with copper-polyethylene shadow cones between the target that 
produces the neutrons and the detector. Figure 2c shows a light discrepancy between 
the measurements and the reconstructed data for the values between the 4 and the 7 
inches spheres; this is probably due to the various efficiencies of the shadowing cones 
for small and big spheres. Anyway the global discrepancy between these two curves 
is 1.2 %. Anywhere else, for the right and left wings, the reconstructed distribution 
is in good agreement with the original data. The unfolded spectrum (fig 2a and 2b), 
clearly shows a mono-energetic spectrum with a mean peaked at 686 keV. The difference 
between 565 keV and 686 keV can be partly explained by the large binning used for 
these energies because of the logarithmic progression of the bins in order to cover a 
wide energy spectrum of nine decades. Also the low energy resolution originated in the 
response matrix can be invoked to explain this discrepancy and finally too a systematic 
bias of about 40 keV from the energy calibration of this Van de Graaf generator. The 
shape of this spectrum is obviously peaked even if the FWHM is quite large : 726 keV. 
The global statistics errors represents 19. % of the unfolding spectrum (see figure 2a) 
while the systematics errors cover 28.8 % (see figure 2b); even if these two figures are 
quite large, this kind of spectrum estimation is sufficient enough for radioprotection 
purposes. These large values of uncertainties are probably due to the low counting rates 
compared to the thermal neutron measurement and also to the bad signal to background 
ratio for direct and diffused neutrons measures. Figure 2d illustrates the high speed of 
the unfolding process: after only 97 iterations, the stability is obtained. 

3.3. Evolution of the minimization criteria Q\ for thermal neutron unfolding 

The progress of the minimization criteria Qx for the unfolding of the thermal neutron 
from Sigma is illustrated in figure 3. The Qx vs. A function is drawn after 1, 2, 3 and 
finally 78 iterations. The evolution of this minimizing criteria is very fast. After the first 
iteration a wide minimum can be seen, while after the second iteration it is already very 
faint and extremly closed to 1. With the evolution of the iteration number the value of 
this minimum get smaller and smaller and also closer to 1. The value of Qi is always 
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+00 this explains why the A minimum is very closed to 1 but never equals 1. So for 
the last iteration before stability, i=78, wc find 1 — A ^ 10~^^ . This kind of unfolding 
scenario is typical of the cases we have studied, and also the order of magnitude of the 
last value of 1 - A lO'^^ 

4. Modification of the interpolation method for the thermal neutron part 
of the spectrum 

4.1. Precise study of the thermal part of the unfolded spectrum of neutrons from the 
Sigma installation 

If we consider into details the unfolded spectrum from the Sigma installation, we 
can be puzzled by the anomalous behaviour of this spectrum in the range 10~^ — 
10~'''Me\^(0.01 — 0.1e\^). Normally, the curve should look hke a MaxweUian distribution 
with a temperature of 300K. 

^ cx E^/^e^kT = 300X = 0.025eV (40) 
dE 

Instead of this behaviour, the neutron spectrum distribution seems to rise very strightly 
at low energies without exhibiting a maximum around 25 meV. Three explanations can 
be argued for this fact. First of all, the low energy cut of the matrix response and 
of the unfolded spectrum is put at 10 meV {10~^eV) which should distort a httle bit 

this low energy region behaviour of thermal spectra. The main reason of this artificial 
cut at 10 mcV is due to the very bad uncertainties of the response functions of the 
Bonner Sphere system to very low energy neutrons produced by the MCNP montc- 
carlo program in the range ImeV - 10 mcV. So, very large statistical error bars are 
obtained through this method at 1 meV (bigger than 20%). This bad quality of results 
produced by MCNP for low energies neutrons is originated by the importance of the 
kinetic effect that polyethylen for example can generate through thermal exchanges 
between the molecules of polyethylene, whose agitation and binding energy are greater 
than 22 meV, and the impinging neutrons of very low energy 1 meV. Thus, these very 
cold neutrons are heated by the brownian motion of the molecules of the detector and a 
good deal of these neutrons do not penetrate inside the detector but are preferentially 
bouncing at the first contact with the surface of the sphere. This explanation is also at 
the origin of the discontinuities observed on the matrix response while the diameter of 
the spheres varies between the bare detector and the first 3 inches sphere for energies 
lower than 1 eV (see figure 4). The bouncing neutrons, due to thermal exchange between 
them and polyethylene explains the reheating and freezing neutrons between 1 meV and 
1 eV. It produces for the spheres of diameter ranging between 1.3 to 2.5 inches at some 
energies an enhancement of the detector's efficiency and a reduction at other energies 
(see figure 5). This last point brings serious consequences while big discontinuities 
appears at low energies for small spheres between naked counter and the 3 inch sphere: 
we are not allowed to perform an interpolation of the measured data for hypothetical 
spheres in the range of the bare detector up to the 3 inch sphere. (Note : diameter of 
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the naked counter =1.29 inch) This is equivalent to say that the measurement function 
is not smooth enough for spheres with diameters lower than 3 inches to authorize a 
direct interpolation. A solution would be to perform some more measurements with 
spheres of lower diameters than 3 inches, but due to these heavy discontinuities, small 
variations in diameter (less than 1mm) on the polyethylene thickness should generate 
strong variations on the neutron's rates. 

4-2. Rectifying the interpolation 

We propose an alternative solution to this problem : in most of the actual cases, the 
low energy spectrum follows a Maxwellian distribution; the very rare cases where this 
assumption is wrong are very few and are well known situations as for example : a specific 
low energy neutrons filtered beam produced with a nuclear reactor adapted to solid 
state physics and cristallography. In other words, guessing a Maxwellian distribution is 
a low cost, extremly realistic and highly efficient assumption for low energies neutrons 
encoutered in actual radioprotection situations. So we performed an interpolation of the 
measurements for the spheres between 2.9 inches up to 12 inches, while we generated 
for a pure maxwellian distribution response of the spheres lower than 1.9 inches down to 
the naked counter, matching the data measured with this bare counter and its response 
to this pure MaxweUian distribution. For the hypothetical spheres between 1.9 and 2.9 
inches, we performed a linear transition to connect these two curves of reconstructed 
measurements versus diameters of the spheres. Unfolding with this method the Sigma 
spectrum shows a satisfying solution regards to the Maxwellian distribution. The price 
to pay is a loss of continuity between the thermal and the epithermal parts of the 
unfolded spectrum. Anyway, the epithermal part is rather small compared to the pure 
thermal part and secondly the incidence of this epithermal part is rather small on the 
total dose flux measurement, because it deals with few neutrons ranging between 1 and 
100 eV. The interesting point of this method is that it doesn't affect at all the unfoldings 
of greater energies neutrons spectra, for example the 565 keV neutrons monoenergetic 
spectrum is absolutely unchanged! 

4-3. Revisiting the Sigma thermal neutron unfolded spectrum in the light of the 
interpolation rectification 

Using the rectification of the interpolation method with a Maxwellian distribution 
forcing for the thermal part of the spectra, onto a real data set from the Sigma 
thermal neutron facility, Am-Be sources + graphite moderator (see figure 6c). We 
obtain a very interesting neutron unfolded spectrum which exhibits a clear Maxwellian 
distribution for thermal neutrons of energy below 10~^e\/(100meK), a high energy 
part of the spectrum for neutrons ranging from IMcV up to 10 McV corresponding 
to direct neutrons produced in the Am-Be source (see figure 6a). Between 10 and 100 
eV (10~® — lO^^MeV"), we can find an enhanced epithermal band, that seems to be 
disconnected from the Maxwellian thermal distribution. Here is the price to pay for this 
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linear combination of two method of data interpolation (Spline+Maxwellian forcing) for 
the small diameters sphere countings. Anyway, this small epithcrmal band is not of a very 
importance for the dose flux calculation. Also, the global statistics errors remains very 
low: 0.24% , compared to 0.11% previously(see figure 6a). While the global systematics 
errors is slightly increased from 3.6% to 5.15% (see figure 6b). We can see that the final 
distribution is obtained after 85 iterations with a minimizing parameter 1 — A less than 
10~^^. The modification of the interpolation gives more realistic results when the original 
approach shows a small deficience due to discontinuities while interpolated. Fortunately 
we have a very good knowledge, thanks to the thermodynamic of the standard thermal 
distribution of low energies neutrons. 

5. Conclusion 

The development of this unfolding iterative Tykhonov approach was leaded by three 
important requirements : a high efficiency and fast converging method, a better and 
faster convergency than the least square method or than the classical Tykhonov method, 
but with the stability of the solution toward fluctuations of the measurements produced 
by the Tykhonov approach and also, a complete unsensitivity of the solution with respect 
to the initial guess. The iterated Tykhonov unfolding approach is modiflcd to ensure a 
good convergence with a parameter running between and 1. An accelerating algorithm 
is introduced between two classical Tykhonov iterations, this algorithm looks like a 
gradient method with a least square term. The calculation of statistics and systematics 
uncertainties is studied in details. The efficiency of this approach is presented in 
two experimental neutron spectrum unfoldings. This method is fully dedicated to 
Bonner Sphere system unfolding. Due to its very wide energy coverage, the Bonner 
sphere system is a very usefuU device for radioprotection evaluation, but it cannot 
be used as a high resolution neutron spectrum instrument for metrology, even with 
an excellent unfolding method! One of the very big advantages of the method here 
illustrated is that no a priori information is introduced during the unfolding process 
leading to coherent results with the expected spectra. This independancc toward the 
initial guess is extremely important for radioprotection measurement when a-priori no 
information on the neutron spectrum is available; not observing this point should lead 
to very prejudiciable wrong evaluations of dose flux calculation which is not allowed 
for radioprotection expertise approach. Nevertheless a very slight stretch to this law is 
applied in such a way that the interpolated counting rates between naked counter and 
the first 3 inches sphere is assumed to correspond to a pure Maxwcllian distribution 
for the thermal neutrons part. This rectification is done to prevent from the high 
discontinuities of the matrix response for low diameters and low energies. This type of 
unfolding gives realistic results for thermal neutrons without modifying the high energies 
neutron spectra. 
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Figure 1. Sigma spectrum unfolding 
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Figure 2. 565 keV spectrum unfolding 




Figure 3. Evolution of Qx vs. A with iteration number 
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Figure 4. Matrix Response of Bonner spheres versus logio{NeutronEnergy /Mev) 
versus Sphere diameter (inch) 
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Figure 6. Sigma spectrum unfolding with modified interpolation : linear combination 
of Spline + Maxwellian distribution forcing in the thermal part 



